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Electronic eigen-states of a square graphene quantum dot(GQD) terminated by both zigzag and 
armchair edges are derived in the theoretical framework of Dirac equation. We find that the Dirac 
equation can determine the eigen-energy spectrum of a GQD with high accuracy even if its size is 
reduced to a few nanometers. More importantly, from the Dirac equation description we can readily 
work out the number and energy gap of the conjugate surface states, which are intimately associated 
with the magnetic properties of the GQD. By using the Hartree-Fock mean field approach, we study 
the size dependence of the magnetic ordering formation in this square GQD. We find that there 
exists a critical size of the width between the two zigzag edges to indicate the onset of the stable 
magnetic ordering. On the other hand, when such a width increases further, the magnetic ground 
state energy of a charge neutral GQD tends to a saturated value. These results coincide with the 
previous results obtained from the first principle calculation. Then, based on the Dirac equation 
solution about the surface state, we establish a simple two-state model which can quantitatively 
explain the size dependence of the magnetic ordering in the square GQD. 

PACS numbers: 75.75.+a, 73.20.-r, 75.50.Xx, 75.70.Cn 



I. INTRODUCTION 

Graphene has become a subject of intense interest since the experimental success in fabricating such 
an atomically thin layer of graphite The valence electron dynamics in such a truly two-dimensional 
material is governed by a massless Dirac equation. As a result, graphene exhibits many unique electronic 
properties03> comparison with the conventional semiconductor materials. From an application point of 
view, graphene possesses very high mobility even at room temperatureQ. Moreover, the planar geometry 
of graphene is of advantage to tailor various nanostructures by the current experimental means, such as the 
lithographic techniques^. So far, the obtainable graphene nanostructures include the one-dimensional(lD) 
nanoribbonsfG^ and zero-dimensional quantum dots[T|. These nanostructures are viewed as the elemental 
blocks to construct the graphene-based nano-devices. 

Accompanying the extensive investigations on the electronic properties in bulk graphene. Graphene nanos- 
tructures also draw much attention of theoretical study P. First of all, some theoretical approaches to 
produce the effective electron confinement in graphene were proposed [§, [13, [111, [l^l which is a nontrivial 
problem due to the Klein tunneling of the carrier in graphene[13l Il4|. For example, GQD structures can be 
formed by patterning gates on a semiconducting graphene nanoribbon[Tol [Tl|. or by using inhomogeneous 
magnetic fields Then, some device applications of the graphene nanostructures were suggested, such 
as the spin qubits based on the coupled GQDs[l0'|. In addition, some electronic properties of graphene 
nanostructures are expected to be different from bulk graphene because of the quantum confinement and 
the edge effect. For instance, the spontaneous magnetization is anticipated to emerge in some graphene 
nanostructures (isl Il6j . which is attributed to the spin polarized electron occupancy at the zigzag- type edges 
of the nanostructures. And such a magnetic ordering has been experimentally demonstrated (ITi, il8|. Quite 
recently, the possible i nag netism of graphene nanostructures with different shapes is theoretically studied 
in some detailsfiol [20l l2ll. [23, [2^ . For example, an infinitely long graphene nanoribbon with zigzag edges 
is possible to behave as a half-metallic material, in which a spin polarized current can be formed[25|. Such 
a property is controlled by an external electric field, which can tune the asymmetry of the band structures 
of the opposite spin electrons. Apart from the ID graphene nanoribbon, GQDs with different shapes also 
exhibit magnetic orderings, such as sqiiare quantum dot and triangular and hexagonal quantum dots ter- 
minated by zigzag edges|19l. [20l [2ll [23 |. Theoretical calculations indicate that in these zero-dimensional 
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graphene structures the magnetic ordering is so robust that can be detected at room temperature. In addi- 
tion, for the square and hexagonal (2Ck i21j quantum dots, there is a critical size which marks the onset of the 
spin-polarized ground state. On the contrary, for a quantum dot with a size smaller than the critical size, 
the ground state is a paramagnetic state. 

The aforementioned results about the magnetic properties of graphene nanostructures are obtained by 
means of the first-principle calculation [2l| as well as the mean field approximation on a Hubbard model of 
a hexagonal lattice. Although these theoretical approaches can give the reliable results about the electronic 
properties of graphene nanostructures, their applicability is restricted within those structures with relatively 
small size. Furthermore, in most cases these approaches can not provide us with a clear physical picture to 
explain the numerical results. For example, an unambiguous explanation about the critical size for the onset 
of the magnetic ordering in GQDs is yet lacking. On the other hand, the Dirac equation description can 
just compensate for the disadvantage of the two theoretical approaches mentioned above. As a theoretical 
model based on the effective mass approximation, a massless Dirac equation can well describe the electron 
properties of the bulk graphene as well as the graphene nanostructures with relatively large sizes [13, HI]- 
Usually, such a model can afford analytical results, which are very helpful to explain intuitively the electronic 
properties associated with the relativistic quantum mechanical feature of graphene. 

So far, the Dirac equation succeeds in describing the band structures of graphene nanoribbons with distinct 
edge types [U, In the present work, we will employ this model to study the electron states in a square 
GQD. With an appropriate boundary condition, we can derive an analytical solution of the electron eigen- 
state in such a GQD. Moreover, by a numerical calculation carried out from the mean field approximation 
of the Hubbard model, we investigate the size dependence of the magnetic property of the GQD. We find 
that there are not only the critical size, but also another characteristic size to indicate the saturation of the 
magnetization. Namely, it is a relatively larger size than the critical size, beyond which the spin polarization 
energy no longer varies with the further increase of the size of the GQD. Then, based on the analytical 
result obtained from the Dirac equation, we establish a simple theoretical model, which can not only reveal 
the physical nature of the emergence of the critical and saturated sizes, but also provide an simple way to 
rapidly create the quantitative result about the critical and saturated sizes, well agreeing with the numerical 
result of the mean field approximation. 

The rest of this paper is organized as follows: Starting from the Dirac equation and using appropriate 
boundary conditions, the wavefunction and the dispersion relation of the electron eigenstate of the square 
GQD are derived in section |TT1 Then the eigen-cnergy spectra calculated from the Dirac equation and the 
tight binding model are compared. In section IIIIl the magnetic property of the GQD is investigated by 
means of Hartree-Fock mean filed theory. By establishing a two-state model, the size dependence of the 
magnetic ordering is quantitatively explained. Finally, the main conclusion is briefly summarized in section 



II. THE ELECTRONIC STATES OF GQD 



The honeycomb lattice of the square quantum dot made of graphene monolayer is schematically shown in 
Fig[TJ The edges of the square GQD are of two kinds, zigzag edges at the top and bottom, and armchair 
ones at the left and right sides. We assume that the dangling a bonds at the edges are passivated by 
hydrogen atoms. Thus, the behavior of the tt band electron near Fermi level is not nontrivially affected 
by the truncation of the a bond at the GQD edges. Within the effective mass approximation, the envelop 
function of the tt band electron in graphene monolayer obeys the following Dirac- like equation [26| 
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where j=\/3tao/2 with t being the nearest neighbor hopping energy. In what follows we use units such that 
7 = ft = 1. k±=kx ± iky and /^^.(j,) = —idx{y) is an operator to measure the momentum deviation from 
K — (— 47r/3ao,0) or K.' = (47r/3ao,0) point. The four components of the spinor wavefunction in Eq.([T]) 
are associated with the total wavefunction ^'(r) by the following relationship. 
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FIG. 1: A schematic of the honeycomb lattice of a square GQD terminated by both armchair and zigzag edges. There 
are two type of atoms(» :A and o : B) The realistic boundary of the GQD is denoted by the thick lines. A hard wall 
labeled by the peripheral rectangular framework consists of the carbon atoms just disconnected from the atoms at 
the edges of the GQD. The electron probability amplitudes at the atoms on the Hard wall must vanish. The size of 
GQD in X direction is denoted by N ox L'^ = {N + l)ao/2, and in y direction by Af or L'y\^ao{M/2 + 1/3). ao is 
the lattice constant. The upward and downward arrows at the opposite zigzag edges denote the net spin moments, 
which indicates the anti-ferromagnetic state. On the contrary, if the net spin moments at the two zigzag edges point 
at the same direction, it indicates the ferromagnetic state. 



where S,(r — JS^) is the carbon atomic wavefunction centered at i?^. And ^fi{R^) denotes the probability 
amplitude of the valence electron appearing in the vicinity of this carbon atom. For the bulk graphene, by 
solving Eq.([T]) we can obtain the electron eigenstate which has the linear dispersion relation e — E/^ — sk, 
and s=±l denoting the conduction and valence bands, respectively. 

As for the present square GQD structure, the tt band electron obeys the same Dirac equation, but is 
subject to the following boundary conditions. At the zigzag edges 



0B(y = 0) = = 0) = q^Aiy - K) = <P'Ay = L') = 0, 



and at the armchair edges 



0) - c^'Jx = 0), - L',) = (x = L',). 



(4) 



(5) 



These boundary conditions have been successfully used to work out the band structures of the graphene 
nanoribbons with different edges [2^ and the energy spectrum of GQDsfl^, [llj. They originate from the 
requirement that the electron probability amplitude at the hard walls around GQD must vanish. Combining 
Eq.([T]) with Eqs.(l4][5|), we can derive the eigen solution of electron states in the square GQD. It is given by 
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The corresponding eigen-energy is given by 



En = S\/pl + q^- 



(6) 



(7) 



Although it takes the same form as the dispersion relation of the bulk graphene, in the present square GQD 
the wavevectors p„ and q in x and y directions are both discrete. p„ is given by 
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Corresponding to a given p„, q is determined by a transcendental equation 

q =p„tangL'y. 
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By analyzing the above equation we find that the electron state with an imaginary wavevector 5 = is 
allowed in the region p„ > 1/-^^, In such a state the electron wavefunction is localized in the vicinity of the 
zigzag edges. Accordingly, it is named as a surface state. In contrast to the surface state, we call a state 
with a real q as the confined state. The normalization coefRcient of the spinor wavefunction for the confined 
state is 

C = (2L;L^ - sin (2p„L^)L;/p„)-i/2, (10) 

whereas for the surface state it is given by 

C = ((2L^ - sin (2p„L^)/p„)(smh {2qL'y)/q - 2L;)/2)-i/2. (n) 

The eigenstate of the GQD shown in Eq.® can be understood in the following way. If the wavefunction 
is rewritten in a form $ = <i>^ + <I>j^?" with $^ = ^b, 0, 0]^ and <^~^" = [0, 0, </)^, 0^]^, we can 
immediately find that and <i']^?" are just the wavefunctions of the eigenstates of an infinitely long zigzag 
nanoribbonjl^] with the free wavevectors p„ and —pn respectively. Therefore, the eigenstate of the GQD 
consists of the linear combination of the eigenstates of the zigzag nanoribbon in K and K' valleys. The 
boundary condition in the armchair direction admixes K and K' valleys. According to such an argument, 
we should restrict the possible surface states in the wavevector range l/i^ < Pn 1^ ""/Sao, where 7r/3ao is 
just the midpoint between K and K' valleys in x direction. When p„ is beyond such a value, there is no 
longer any new surface state due to the valley admixing. Considering such a wavevector limit, we can readily 
determine the number of the surface states by simply counting the number of the discrete wavevectors p„ in 
this range. For example, for a GQD of size TV = 13 and M = 10(denoted for short as N13M10), the allowed 
Pn are 7r/21ao and 47r/21ao. Therefore, the total number of the surface states in valence and conduction 
bands are 2Nt = 4, where Nt is the number of the allowed p„ in the range of < Pn < 7r/3ao. Moreover, 
from this method we can infer that when TV < 7 no surface state survives, independent of the value of M . 
We will see below that the surface states are responsible for the magnetic property of the GQD. Such a 
simple way to determine the number and energy of the surface states is helpful for us to explain intuitively 
the numerical result about the magnetic property of the GQD. In addition, it should be noted that in the 
following discussion, we only consider the GQD with odd N and even M. Such a case corresponds to a GQD 
without any carbon atom at the boundary connected to the GQD by a single tt bond. Finally, although the 
Dirac equation can give an analytical description of the electron eigenstate in the square GQD, we have to 
point out that its applicability should be strictly restricted within the linear dispersion region of the tt band 
of graphene. It is known that the energy scope of the linear dispersion region in the n band is about i/3, 
away from the Dirac point. Accordingly, in the two dimensional k-space the linear dispersion region forms a 
circle around the K or K' point with a radius equal to 2/{3^/3a). Therefore, if the energy of an eigenstate of 
the square GQD is much lower than t/3, it can be well described by the Dirac equation. In contrast, when 
an eigen-energy exceeds i/3, the Dirac equation gets poorer to describe such an eigenstate in the square 
GQD. From Eq.® we can see that the interval between two adjacent wavevectors in the armchair direction 
is Apn = tt/L^. The discrete wavevectors in the zigzag direction should be numerically determined from 
Eq.([n|). Therefore, it is not straightforward to find the interval between the adjacent wavevectors in this 
direction. However, for a relatively small Pn we infer from Eq.([9]) that such a wavevector interval is roughly 
equal to Aq = / ^'y Thus, the two kinds of wavevector interval are inversely proportional to the sizes in 
the respective directions. When Ap„ and Aq become comparable to the diameter of the circle of the linear 
dispersion limit, there is hardly any eigenstates in the linear dispersion region. Thus, from the relation 
Ap„ ^ Aq^ 4/(3A/3a) we can find the minimal size of the GQD for the complete invalidity of the Dirac 
equation description. By a simple evaluation, we find that such a minimal size is Nmm = 7 and Mmm = 8. 
This means that the GQD with size of N7M8 is the smallest one to which the Dirac equation is applicable. 

To check the validity of the Dirac equation solution about the electronic eigenstate of the GQD in 
some details, we compare the low- lying energy levels calculated by solving Eqs.([7][5]) to the ones ob- 
tained from the tight-binding model. The tight-binding Hamiltonian of the square GQD takes a form 
as Hfb = —tJ2{ij) J^ai'^iAa'^jBo- + h.c) where c\^^ is the electron creation operator associated with a local 
atomic state at lattice point i; And denotes any pair of the nearest neighboring carbon atoms, a —] (|) 
corresponds to the up and down spins. In the basis set consisting of the local atomic orbits, the tight-binding 
Hamiltonian changes into a matrix. By diagonalizing this Hamiltonian matrix we can obtain the electronic 
eigen-energy spectrum and the eigen wavefunctions. Noting that electronic eigen-states are spin-degenerate, 
though we write the spin index explicitly in the above tight-binding Hamiltonian. The electron spin becomes 
relevant only in the self-consistent calculation of the electron energy spectrum in the next section where the 
Hubbard interaction is taken into account. 

The comparison of the numerical results of the low-lying eigen-energy spectra obtained by the Dirac 
equation as well as the tight-binding model is visualized in Figl^] For a relatively large GQD, N35M24 as 
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FIG. 2: A comparison of the low-lying energy levels calculated from the Dirac equation (triangle,cross) and tight- 
binding model (square) for GQDs with different sizes, (a) N35M24; (b) NI1M24; (c) N7M24;(d) N7MI2; (e) N7M8 
and (f) N5M4. 
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FIG. 3: Low-lying energy levels versus pn for a GQD of size N13M10 calculated from Dirac equation (triangle). 
The dashed lines depict the dispersion relation of an infinitely long zigzag nanoribbon with width M=10, which is 
calculated from the tight-binding model. The inset shows more clearly the two pairs of the conjugate surface states of 
the GQD. And the energy gaps relative to the largest p„ and the smallest p„ are denoted by Eg and Eg respectively 



shown in Fig|21(a), the Dirac equation result agrees with the tight-binding result very well. Then we reduce 
the size of the GQD only in x direction. These results are plotted in Fig[2lja)-(c). We can see that the Dirac 
equation results get poorer with the decrease of N. In Figl^Jc), despite N = 7, there are a few low-lying 
eigen- energies obtained by the Dirac equation solution to be close to the tight-binding result. This is due 
to that the relatively large M retains these eigenstates in the linear dispersion region. When M decreases 
further, the results calculated by the two models deviate from each other notably, as shown in Figl^fc)- 
(e). The result shown in Figl^I^f) demonstrates that the Dirac equation description fails to give the correct 
eigen-energy spectrum for the GQD smaller than the one of N7M8. On the other hand, the result shown 
in Fig[2fa) indicates that at least 10 low-lying eigenstates can be safely described by Dirac equation for a 
GQD with size of 4nmx5nm. In conclusion, the numerical comparison made in Fig|2] supports our simple 
criterion given above for the applicable limit of the Dirac equation approach to the GQD. 

In Fig[3]some low-lying eigen-energies versus p„ are plotted for a GQD of size N13M10, which is compared 
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with the dispersion relation of a zigzag ribbon with width M=10. Although these energy levels of the GQD 
are discrete, from this figure we can readily infer that the energy-wavevector relation of a GQD will change 
into the band structure of a zigzag nanoribbon with the continued increase of the size N. Thus, we can say 
that the dispersion relation of the zigzag nanoribbon remains in the Dirac equation description of the square 
GQD. This is one advantage of the Dirac equation over the tight-binding model in describing the electron 
states of the GQD. In addition, from the Dirac equation solution, we can easily distinguish the surface 
states from the confined states, because the two kinds of state have distinct forms of wavefunction. This can 
be viewed as another advantage of the Dirac equation description. In contrast, it is difficult to identify a 
surface state in the tight binding model, in particular, for a GQD with small size. In fact, the Dirac equation 
method was previously used to describe the electron states in other GQDsfl^, El, For example, for a 
GQD formed by applying a gate voltage on a graphene nanoribbon the surface states and the confined 
states can also be easily distinguished from each other by using the Dirac equation method. Finally, from the 
inset of Fig[3]we can clearly see that corresponding to any allowed p„, there are two conjugate surface states, 
belonging to the valence and conduction bands respectively. And between them there is a finite energy gap. 
Although there are only two pairs of surface states for the GQD shown in FigEl a common feature about 
the surface state visualized in this figure is that the pair of surface state with the minimal Pn has just the 
maximal energy gap and vice versa. 



III. MAGNETIC PROPERTIES OF THE GQD 



Magnetic properties of graphene nanostructures with various shapes have drawn considerable interest. For 
example, some theoretical investigations indicate that a zigzag nanoribbon has a spin-polarized ground state, 
which is tightly associated with the surface state at the zigzag edges. This implies that a spontaneous mag - 
netization may occur in such a nonferromagnetic material. Motivated by these previous workfisl [20l. [2ll. [2^ . 
we now study the possible magnetic property of the square GQD. To do this, we adopt a single band Hubbard 
model(to incorporate the Hubbard terms into the tight binding model), and treat it within the Hartree-Fock 
approximation. It was previously proved that most magnetic properties of graphene nanostructure can be 
captured by such a simple approach[l^ [l^, The Hartree-Fock mean field Hamiltonian of the present 
GQD take a form 

H = Htb + U ^(("vi)ni^T + ("■vT)"-vi ~ ("vT>("-vi>)- (12) 

where ?T.i^|(= cl^t|Ci/jT) ^^'^ ("-vt) electron number operator and the average electron occupation at 

an arbitrary lattice point respectively. Besides, U is the on-site Hubbard energy. For a charge neutral GQD, 
the single-electron picture adopted above tells us that all eigenstates in the valence band are fully occupied, 
whereas all states belonging to the conduction band are empty. However, the finite Hubbard U distorts such 
a simple electron distribution since it resists the double occupancy of a surface state by two opposite-spin 
electrons. As a result, spin polarized electron occupancy on individual lattice points may occur in the charge 
neutral GQD. 

The eigen solution of the Hartree-Fock Hamiltonian for a charge neutral GQD can be obtained by iteration 
method. In analogy with the previous work, our iterative calculation indicates that the spin polarization 
situation of the obtained eigen-state of the charge neutral GQD depends on the initial spin configuration 
to start the iteration procedure. The different initial states will lead to the eigen-states with distinct kinds 
of spin polarization. At first, to begin with an initial spin configuration of Neel order, we will arrive at an 
eigen-state with spin polarized electron occupation on individual lattice points. In particular, on the lattice 
points near the two zigzag edges the spin polarization is very strong. The net spin distributions at the two 
zigzag edges show the anti-ferromagnetic order(AFM) in FiglT] Second, if we start from an initial state 
with the uniform spin polarization at all the lattice points, the self-consistent calculation converges to an 
eigen-state with ferromagnetic(FM) ordering. In addition, a paramagnetic state(PM) can be achieved if the 
initial spin configuration is set to be unpolarized at all lattice points. Herein we adopt the same definition 



about the magnetic orderings as given in the previous works [15|. The AFM state refers to that the spin 
moments of the carbon atoms on one zigzag edge are anti-aligned to that on the opposite edge, while the 
FM state means that the spin moments on both zigzag edges point at the same direction. Moreover, the 
PM state can be defined alike, which refers to that the spin up and down electrons are equally occupied at 
every lattice point. Next we focus on the size dependence of the magnetic orderings. To do this, in Figd] 
we compare the total energy difference(Ai5) between the FM(AFM) states and the PM state for the charge 
neutral GQD, as a function of the width M between the two zigzag edges, while the width N in the armchair 
direction takes several typical values. First of all, we can find that there exists a critical size Mc, which 
denotes the onset of the PM-FM(AFM) transition. Namely, only when M > Mc is the magnetic ordering 
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FIG. 4: Energy difFerences versus M between AFM(circle) and PM as well as FM (triangle) and PM states for charge 
neutral GQDs with diflterent sizes and on-site energies. (a)N7U0.1t; (b)N7U0.5t; (c)NllUO.lt; and (d) NllUO.St. 
The critical size Mc and the saturated size Ms are labeled on the curves. 



stable. By comparing the results shown in FiglDJa-d), we find that Mc depends on the transverse width iV 
and the Hubbard U sensitively. With the increase of iV and C/, the critical size Mc decreases notably. 

There is one point we have to emphasize herein about the magnetic ordering in the GQD depicted above. 
The FM or AFM state in the GQD considered by us and the one in an extended honeycomb lattice have 
distinct underlying mechanisms. It is the surface states localized at the zigzag edges to cause the AFM in 
the GQD. The surface states in the GQD tend to dispersionless as the width between the opposite zigzag 
edges getting large. Thus, a finite density of states forms in the vicinity of the Dirac point(energy zero 
point), which is responsible for the formation of the magnetic ordering even with a very small Hubbard U. 
However, such a surface state is absent in the extended honeycomb lattice in which there is no zigzag edge. 
The linear dispersion of the extended honeycomb lattice leads to the vanishing density of states at the Dirac 
point, which then requires a very large Hubbard U for the PM-AFM transition in such an extended lattice. 
In fact, a similar comparison was made in Refd where the magnetic ordering in a zigzag ribbon only needs 
a much smaller (infinitely small, in fact) Hubbard U than in the extended honeycomb lattice. In analogy 
with the GQD, the zigzag ribbon possesses surface states which form a flat band at the Dirac point, leading 
to the magnetic ordering even with an infinitely small U. In short, we think the magnetic ordering formed 
by a small U is due to the existence of the surface state in the GQD. The spin polarized electron occupancy 
is notable only in the regions near the zigzag edges. The value of the Hubbard U we used in our work is 
too small to show the formation of the AFM state following the mechanism of the extended lattice(from the 
result shown in Ref[li we can find out that such a critical U is larger than 2t). 

Another noticeable feature in Figdjis that the energy difference tends to a saturated value when the width 
M exceeds a specific value M^. Hereafter we call M^ the saturated size. It is quantitatively determined 
following such a way: when the size of the GQD increases from Mg to Mg + 2, the relative increment of A_E 
should be less than 1%. In addition, from Fig|4]we can find that the energy of the AFM state is notably lower 
than that of the FM state in the region Mc < M < Mg. This indicates that the AFM state in the region 
between the critical size and the saturated size is just the ground state of the charge neutral GQD ^.1271. 
On the other hand, when M is sufiiciently large, the energy difference between the FM and AFM states 
becomes indistinguishable. This can be explained in such a way: The surface states in the GQD can be 
viewed as the bonding or anti-bonding states arising from the interaction between two kinds of surface states 
located at the opposite zigzag edges. In fact, each of such two surface states belongs to a semi-infinite 
two-dimensional graphene terminated by a zigzag edge. When M is very large, the interaction between the 
two kinds of surface states gets weak. As a result, the exchange integral between them which determines the 
relative orientations of the net spin moments at two zigzag edges becomes negligibly small. Thus the energy 
difference between FM and AMF states tends to zero. Finally, we have to point out that the existence of the 
critical size for a hexagonal and square graphene quantum dots was previously reported [1^ [2l], (2^, based 
on the first principle calculation. The result about the square GQD agrees quantitatively with our present 
result obtained by the mean field approximation (2ll. [2^ . However, a clear explanation about the critical size 
is yet lacking. 

Next we try to give a reasonable explanation about the occurrence of the critical and saturated sizes. 
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At the first step, we discuss the relation between the magnetic order formation and the electron occupancy 
on the surface states. Accordingly, we can work out a simple criterion which can qualitatively explain the 
numerical result about the critical size obtained above by the Hartree-Fock mean field theory. Then, we 
establish a two-state model which can give a quantitative explanation to the numerical result about the 
critical and saturated sizes. All these arguments benefit from the Dirac equation solution about the single- 
particle surface states. As analyzed in the preceding section, the number of the surface states in the range 
< Pn < 7r/3ao is finite. Corresponding to a specific p„, there are two conjugate surface states, belonging 
to the conduction and valence bands respectively. And the energy gap between them depends on p„ and the 
size of the GQD. We consider specially the two conjugate surface states with the maximal p„. In comparison 
with other surface states, this pair of surface states has the smallest energy gap. The finite Hubbard U alters 
the single electron energy spectrum since it affords an on-site Coulomb repulsive potential. For example, 
in a presumed paramagnetic state, the surface state with a specific spin a in the valence band will rise by 
U{nia) = U/2 due to the Coulomb repulsion. When such a shift makes the surface state in the valence 
band is aligned with the surface state in the conduction band, the system is likely to show the spontaneous 
magnetization. Following such an analysis, we obtain a simple criterion to determine the critical size. Mc is 
the size at which the inequality 

Eg < U/2 (13) 

begins to hold true, where Eg is the energy gap between the pair of surface states corresponding to the 
maximal p„. By calculating Eg we can obtain the critical size M^. Using such a simple criterion, we 
estimate the critical size Mc, varying with the size N as well as the Hubbard U. These results are shown 
in FiglHl In comparison with the Hartree-Fock mean field result, we find that the criterion given by Ea. (jl3p 
can roughly account for the dependence of the critical size on N as well as U. 

Although the above criterion is not too bad, we will show that a more quantitative explanation about the 
critical size is available. Now that the onset of the magnetic ordering is controlled by the pair of surface 
states with the minimal energy gap, we establish a simple two-state model by retaining only this pair of 
surface states in the mean field Hamiltonian. Thus, the Hartree-Fock Hamiltonian for this two-state model 
takes a form as 

H*" = e.cl^c,^ + U Y{m^^,a)mif,^■ (14) 

here s = ± denotes the two surface states belong to the conduction and valence bands respectively. Csa{c\c,) is 
the electron annihilation(creation) operator of the surface state of quantum index sa. The electron number 
operator rhi^rr counts only the contributions of the two surface states to the electron occupancy on the 
individual carbon atoms. According to such a meaning, it is associated with the electron number operator 
of the two surface states via a relation 

m^^,a=Y\'ipl,{^)?^cl„Csa, (15) 

here J7 = is the area of unit cell, and liA^lOP is just the probability of the electron in surface state 

(s) appearing in the vicinity of the carbon atom /i at lattice point i, which can be directly calculated from 
the analytical wavefunction in Eq.Q. For the charge neutral GQD, the two surface states accommodate two 
electrons with opposite spins. We thus have X^so- ^la'^sa = 2. In addition, it should be noted that in such 
a two-state model we have ignored the contributions to the average electron occupancy on the individual 
lattice points from other occupied single electron eigen-states, except for the two surface states retained in 
this model. This is because that other occupied electron states in the valence band only provide a spin- 
unpolarized charge background, which infiuences the opposite spin electrons in the two surface states on an 
equal footing. Substituting Eq. p3|) into Eq. p^ . we obtain the following diagonal Hamiltonian 

H'' = Y.ie, + J2 Uo{cl,^c,-^))clc,,, is' = ±s) (16) 

s,a s' 

with 

Uo^J2^-n [ dr\r^{rr-\r^ir)\\ (17) 

Here Uq can be calculated analytically by the Dirac equation wavefunction given in Eq.Q, or numerically 
by tight-binding model instead. By a simple derivation from the Dirac equation wavefunction we obtain the 
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analytical form about Uq- It is given by 

o /qt-t-/ sinh 4qLy sinh2qLy , 3Ly \ 

u r / sinh 2qLy _ o/, )2 ^ ' 

x\ q y J 

For the PM state, two electrons occupy the valence band surface state(e_) with opposite spins. From the 
diagonal Hamiltonian given above, we can immediately obtain that the energy of the two electrons is equal 
to Epm = 2£_ + 2Uo. On the other hand, for the possible magnetic ordering state, the two electrons occupy 
the two distinct suface states with the same spin. The correspond energy is then i?mo = £- + £+■ The 
critical size for the magnetic ordering to become stable corresponds to Emo < Epm, Namely, 

Eg < 2Uo. (19) 

By means of this criterion we can determine the critical size Mc- It should be noted that although the 
two-state model is established according to the Dirac equation description of the surface states, the two 
parameters Eg and Uq can also be calculated from the tight-binding model. Therefore, the two-state model 
is expected to be still valid even when the conjugate surface states are beyond the linear dispersion region. 
The critical size obtained from this two-state model is shown in FigO In comparison with the mean field 
result, we find that the two-state model can give a quantitative explanation about the critical size, as a 
function of U and N. Besides, in FiglSJa) we also find that the two-state model with the parameters 
evaluated from the Dirac equation can no longer predict the critical size satisfactorily with the increase of 
N. This can be readily understood. In the two-state model, we consider the pair of surface states with 
the maximal p„, which gets away from the center of the valley with the increase of A'^. As a result, the 
Dirac equation becomes poorer to give a quantitative description about the electron probability amplitude. 
Instead of the Dirac equation solution, if we evaluate the parameters Eg and Uo from the tight binding 
model, as shown in FigOthe two-state model always gives a satisfactory result, which demonstrates that the 
two-state model has captured the main mechanism dominating the spin polarization in the GQD. Finally, 
we would like to point out that our numerical result about the critical size shown in Fig[5] coincides with 
those obtained by the first principle calculation or the mean field method in the previous workpll [2^ \2§^ . 
In particular, our calculation indicates that when iV < 7 no surface state exists, hence no magnetic ordering 
occurs. And when N = 7 the critical size is Mc — 8. These quantitative results were also produced in the 
relevant work obtained by the first principle calculation [29] . According to solution about the surface state 
given in the previous section, we can predict that the critical size Mc will tend to zero, as the size N becomes 
sufficiently large. This is because that there must be a dispersionless surface state pair when N becomes 
sufficiently large, regardless of the size M. Thus, The GQD becomes an infinitely long zigzag ribbon which 
always possesses the dispersionless surface states. As a result, the spontaneous magnetization can occur at 
an arbitrarily small M. 

Now we turn to discuss the occurrance of the saturated size Ms , based on the Dirac equation description 
of the surface states. The mean field result of Ms as a function of TV is shown in FigEl It depends on N non- 
monotonously, in contrast to its insensitive dependence on the Hubbard U. For example, when N increases 
from 7 to 11, Ms decreases notably, followed by an abrupt rise at iV = 13. Such a situation recurs when 
N increases further. According to our solution about the surface state, when N increases within a small 
range, such as from = 7 to 11, the number of surface states does not change. The pair of surface states 
with the minimal p„ has the maximal energy gap (we denoted it as Eg in the inset of Fig|3]). According to 
our two-state model, such a gap goes against the formation of spin polarization in the two conjugate surface 
states. And the energy difference between the spin polarized and unpolarized electron occupancies on the 
two surface states is equal to Eg — 2Uq. From Eqs.([39]) we infer that with the increase of M, q goes close 
to pn. As a result. Eg tends to vanish. Meanwhile, as q goes close to p„, Uq tends towards a constant, 
because the wavefunction (EqE]) of surface states near zigzag edges tends towards a constant value. When 
Eg becomes sufficiently small, such a pair of surface states, hence all pairs of surface states, have stable 
contribution to the spin polarization. Thus, the total energy difference shown in Fig[4] tends to a saturated 
value. In comparison to the case of = 7, the GQD with N = 11 has a smaller Eg corresponding to the 
same M. Thus the GQD with TV = 11 

corresponds to a smaller Ms. However, when N increases further, a new surface state comes into being, 
with an larger Eg than the case of = 7. Thus a larger Ms is needed. This justifies the abrupt rise of Ms 
at N = 13 as shown in Fig[51 In such a spirit, we can establish a simple criterion by which we can estimate 
the value of Mg. This is 

Pn-q< V- (20) 

where p„ is the smallest wavevector in the surface state allowed region. 77 is an appropriate small quantity 
to characterize the extent that q approaches to p„. The energy gap Eg is then sufficiently small and Uq 
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FIG. 5: (a) The critical size Afc versus the transverse size A''. (b)The critical size Mc versus U. The square symbol 
denotes the Hartree-Fock mean field result; The star symbol is the result obtained from the simple criterion given 
by Eq.(13). The triangle and cross symbols for the results obtained from the two-state model. The former is the 
result that the parameters Eg and Uo are evaluated from Dirac equation. And the latter is the result with Eg and 
Uo calculated from tight binding model. 




FIG. 6: The saturated size Ms versus the transverse size A'^. The Hartree-Fock mean field results for U=0.5t(square) 
and U=t(star) as well as the estimated result by using Eq.20 (triangle) are compared. These results coincide with 
each other very well. 



becomes a constant. For the numerical calculation, we choose = 10^*/ao. From FiglHlwe can find that 
such a simple rule creates a saturated size Ms which agrees with the Hartree-Fock mean field result very 
well. Finally, we have to point out that the result in Fig[6]only shows the saturated size Ms in a very finite 
range of the size N, because that the self-consistent calculation becomes rather time-consuming as N gets 
larger. But we can predict that as N increases furthermore, the oscillation of the Ms will become weak since 
the discrete wavevector P„ tends to a continuous quantity. Finally the saturated size Ms tends to that of 
the infinitely long zigzag ribbon. 
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IV. CONCLUSION AND REMARKS 



The electronic eigenstates of a square GQD terminated by both zigzag and armchair edges have been 
studied analytically in the theoretical framework of Dirac equation. By comparing with the result of tight 
binding model, we find that the Dirac equation can well describe the electron eigen-states even if the size of 
a GQD reduces to a few nanometers. Moreover, the Dirac equation method has advantages over the tight 
binding model in two aspects. At first, the Dirac equation solution about the electron eigen-states can tell 
us not only the energy levels but also the dispersion relation. Then, from the Dirac equation solution, we 
can readily determine the number of the surface states. In addition, by using the Hartree-Fock mean field 
theory, we have also investigated the magnetic properties of the square GQD. We find that stable magnetic 
ordering states are allowed for a charge neutral GQD with an appropriate size. The magnetic ordering 
depends on the width between two zigzag edges sensitively. Only when the width is larger than a critical size 
is the magnetic ordering stable. On the other hand, when this width becomes sufficiently large, the magnetic 
ordering ground state energy tends towards a saturated value. We find that the critical size is dominated by 
the pair of surface states with the minimal energy gap, while the saturated size is determined by the pair of 
surface states with the maximal energy gap. Based on the Dirac equation description, we establish a simple 
model in which only the two dominated surface states are incorporated. Consequently, this two-state model 
can quantitatively explain the size dependence of the magnetic ordering of the square GQD. Thus, by virtue 
of such a toy model, we can estimate rapidly the characteristic sizes for the formation of magnetic ordering 
in the GQD. 
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